timeseries analysis : normality, q-q plots, distr estimation
Andreas Hadjiprocopis
February, 2018
#!/usr/bin/env Rscriptset.seed(1234)
source('lib/UTIL.R');
source('lib/DATA.R');
source('lib/IO.R');
source('lib/TS.R');
source('lib/MC.R');
source('lib/SEASON.R');
source('lib/MIXTURES.R');
infile='cleaned_data/dat1.eliminateNA.csv'
dat1 <- data.frame(read_data(
filename=infile
))## read_data(): data read from file 'cleaned_data/dat1.eliminateNA.csv'.if( is.null(dat1) ){
cat("call to read_data() has failed for file '",infile,"'.\n", sep='')
quit(status=1)
}
dummy <- remove_columns(inp=dat1, colnames_to_remove=c('id'))
dat1_noid <- dummy[['retained']]
dat1_id <- dummy[['removed']]
dat1_detrended_id <- list(c(dat1_id[['id']][2:length(dat1_id[['id']])]))
names(dat1_detrended_id) <- c('id')
# detrend the data
dat1_detrended <- detrend_dataset(inp=dat1_noid,times=1)Check for normality: do data come from a Gaussian distribution process For the Anderson-Darling and Shapiro-Wilk normality test the p-value refers to the null hypothesis being that data DOES not come from a Gaussian distribution process (i.e. high p-value means ‘a high probability that data comes from Gaussian distribution process’) Data (detrended or not) do not seem to come from a Gaussian distribution process because the p-values of the test are very very small. First let’s test a normally distributed sample just to see what kind of p-values to expect.
library(nortest)
print(ad.test(rnorm(n=1000, mean=5, sd=3)))##
## Anderson-Darling normality test
##
## data: rnorm(n = 1000, mean = 5, sd = 3)
## A = 0.57857, p-value = 0.1323print(ad.test(rnorm(n=1000, mean=5, sd=1)))##
## Anderson-Darling normality test
##
## data: rnorm(n = 1000, mean = 5, sd = 1)
## A = 0.18795, p-value = 0.9025print(shapiro.test(rnorm(n=1000, mean=5, sd=3)))##
## Shapiro-Wilk normality test
##
## data: rnorm(n = 1000, mean = 5, sd = 3)
## W = 0.99866, p-value = 0.6587print(shapiro.test(rnorm(n=1000, mean=5, sd=1)))##
## Shapiro-Wilk normality test
##
## data: rnorm(n = 1000, mean = 5, sd = 1)
## W = 0.99902, p-value = 0.8817it looks like if p-value is not extremely small then we can assume it passes the test of normality.
for(acol in names(dat1_noid)){
cat("original data, col", acol, "\n", sep='')
x<-ad.test(dat1_noid[[acol]])
print(x)
}## original data, colA
##
## Anderson-Darling normality test
##
## data: dat1_noid[[acol]]
## A = 75.256, p-value < 2.2e-16
##
## original data, colB
##
## Anderson-Darling normality test
##
## data: dat1_noid[[acol]]
## A = 1147.2, p-value < 2.2e-16
##
## original data, colC
##
## Anderson-Darling normality test
##
## data: dat1_noid[[acol]]
## A = 2030.8, p-value < 2.2e-16
##
## original data, colD
##
## Anderson-Darling normality test
##
## data: dat1_noid[[acol]]
## A = 138.07, p-value < 2.2e-16
##
## original data, colE
##
## Anderson-Darling normality test
##
## data: dat1_noid[[acol]]
## A = 7884.2, p-value < 2.2e-16
##
## original data, colF
##
## Anderson-Darling normality test
##
## data: dat1_noid[[acol]]
## A = 10508, p-value < 2.2e-16for(acol in names(dat1_detrended)){
cat("detrended data, col", acol, "\n", sep='')
x<-ad.test(dat1_detrended[[acol]])
print(x)
}## detrended data, colA
##
## Anderson-Darling normality test
##
## data: dat1_detrended[[acol]]
## A = 3369.7, p-value < 2.2e-16
##
## detrended data, colB
##
## Anderson-Darling normality test
##
## data: dat1_detrended[[acol]]
## A = 3424, p-value < 2.2e-16
##
## detrended data, colC
##
## Anderson-Darling normality test
##
## data: dat1_detrended[[acol]]
## A = 1491.1, p-value < 2.2e-16
##
## detrended data, colD
##
## Anderson-Darling normality test
##
## data: dat1_detrended[[acol]]
## A = 10395, p-value < 2.2e-16
##
## detrended data, colE
##
## Anderson-Darling normality test
##
## data: dat1_detrended[[acol]]
## A = 11831, p-value < 2.2e-16
##
## detrended data, colF
##
## Anderson-Darling normality test
##
## data: dat1_detrended[[acol]]
## A = 12326, p-value < 2.2e-16Conclusion: data do not pass the Normality test. This means that non-parameteric statstical tools should be used for assessing the data. For example Pearson’s require bivariate normality while Spearman’s does not. print summary statistics for original data:
for(acol in names(dat1_noid)){
x <- dat1_noid[[acol]]
cat("Original data, column ", acol, ", mean=", mean(x), ", sd=", sd(x), "\n", sep='')
}## Original data, column A, mean=49.39303, sd=9.499579
## Original data, column B, mean=80.99835, sd=11.71468
## Original data, column C, mean=5.08511, sd=5.262529
## Original data, column D, mean=29.81647, sd=0.3481571
## Original data, column E, mean=83.66461, sd=159.3282
## Original data, column F, mean=0.3703302, sd=0.8356359print summary statistics for detrended data:
for(acol in names(dat1_detrended)){
x <- dat1_detrended[[acol]]
cat("Detrended data, column ", acol, ", mean=", mean(x), ", sd=", sd(x), "\n", sep='')
}## Detrended data, column A, mean=-1.025213e-05, sd=0.2437584
## Detrended data, column B, mean=0.0002468104, sd=1.291171
## Detrended data, column C, mean=0, sd=2.358147
## Detrended data, column D, mean=-2.018302e-06, sd=0.00484357
## Detrended data, column E, mean=0, sd=57.02373
## Detrended data, column F, mean=-9.7284e-21, sd=0.1546259Boxplot of all of original data’s columns in one page for comparison Note: a boxplot shows the minimum and maximum values in the whiskers The interquartile range (between lower and upper quartiles) is the rectangle’s width. The median is the indicator line inside the rectangle.
boxplot(
dat1_noid,
main=paste0('Original data boxplot', sep='')
)
boxplot of original data per column:
for(acol in names(dat1_noid)){
x <- dat1_noid[[acol]]
boxplot(
x,
main=paste0('Original data, col ', acol, sep='')
)
}
boxplot of all of detrended data’s columns in one page for comparison
boxplot(
dat1_noid,
main=paste0('Detrended data', sep='')
)
boxplot of detrended data per column:
for(acol in names(dat1_noid)){
x <- dat1_detrended[[acol]]
boxplot(
x,
main=paste0('Detrended data, col ', acol, sep='')
)
}
Do a scatter plot of pairs of variables with a regression line added for the original data
library(ggpubr)## Loading required package: ggplot2## Loading required package: magrittradf <- list2dataframe(dat1_noid)
cnames <- names(adf)
Ncnames <- length(cnames)
for(i in 1:Ncnames){
acol1 = cnames[i]
for(j in (i+1):Ncnames){
acol2 = cnames[j]
ggscatter(
adf,
x=acol1,
y=acol2,
add="reg.line",
conf.int=TRUE,
cor.coef=TRUE,
cor.method = "spearman",
title=paste0('Scatter plot of original data, columns ', acol1, ' and ', acol2,sep='')
)
}
}Do a scatter plot of pairs of variables with a regression line added for the detrended data
library(ggpubr)
adf <- list2dataframe(dat1_detrended)
cnames <- names(adf)
Ncnames <- length(cnames)
for(i in 1:Ncnames){
acol1 = cnames[i]
for(j in (i+1):Ncnames){
if( j > Ncnames ){ next }
acol2 = cnames[j]
ggscatter(
adf,
x=acol1,
y=acol2,
add="reg.line",
conf.int=TRUE,
cor.coef=TRUE,
cor.method = "spearman",
title=paste0('Scatter plot of detrended data, columns ', acol1, ' and ', acol2,sep='')
)
}
}Print and plot (Spearman) correlation matrix of original data (pairwise). Correlation tells us whether two variables move together in the same or opposite direction (correlations of +1 and -1) or are their movement is not related at all (correlation of 0). Correlation is related to covariance in that it is essentially covariance normalised over the product of the standard deviations of the two variables.
cnames <- names(dat1_noid)
Ncnames <- length(cnames)
amatr <- matrix(unlist(dat1_noid), ncol=Ncnames, byrow=T)
colnames(amatr) <- cnames
library(corrplot)## corrplot 0.84 loadedcor_matrix <- cor(amatr, method='spearman')
cor_pvalues <- cor.mtest(amatr)the correlation matrix of original columns set:
print("correlation matrix of original data:")## [1] "correlation matrix of original data:"print(cor_matrix)## A B C D E F
## A 1.0000000 0.9953697 0.9926465 0.9903925 0.9882165 0.9860988
## B 0.9953697 1.0000000 0.9952715 0.9925822 0.9904287 0.9881424
## C 0.9926465 0.9952715 1.0000000 0.9955148 0.9928157 0.9904968
## D 0.9903925 0.9925822 0.9955148 1.0000000 0.9952246 0.9925382
## E 0.9882165 0.9904287 0.9928157 0.9952246 1.0000000 0.9953121
## F 0.9860988 0.9881424 0.9904968 0.9925382 0.9953121 1.0000000statistical significance of the correlation matrix
print("and its statistical significance:")## [1] "and its statistical significance:"colnames(cor_pvalues$p) <- cnames
rownames(cor_pvalues$p) <- cnames
print(cor_pvalues$p)## A B C D E F
## A 0 0 0 0 0 0
## B 0 0 0 0 0 0
## C 0 0 0 0 0 0
## D 0 0 0 0 0 0
## E 0 0 0 0 0 0
## F 0 0 0 0 0 0acolor <- colorRampPalette(c("red", "white", "blue"))(20)
corrplot(
cor_matrix,
p.mat=cor_pvalues$p,
method="ellipse",
order="original",
col=acolor,
type='upper',
sig.level=0.0,
insig="p-value",
title='correlation matrix of original data'
)
plot correlation matrix of detrended data (pairwise)
cnames <- names(dat1_detrended)
Ncnames <- length(cnames)
amatr <- matrix(unlist(dat1_detrended), ncol=Ncnames, byrow=T)
colnames(amatr) <- cnames
library(corrplot)
cor_matrix <- cor(amatr, method='spearman')
cor_pvalues <- cor.mtest(amatr)the correlation matrix of detrended columns set:
print("correlation matrix of original data:")## [1] "correlation matrix of original data:"print(cor_matrix)## A B C D E F
## A 1.00000000 0.01846806 0.08630014 0.10118207 0.07398528 0.06682168
## B 0.01846806 1.00000000 0.01899633 0.09033755 0.10640717 0.08220713
## C 0.08630014 0.01899633 1.00000000 0.02234015 0.09594644 0.10275620
## D 0.10118207 0.09033755 0.02234015 1.00000000 0.02221574 0.08849447
## E 0.07398528 0.10640717 0.09594644 0.02221574 1.00000000 0.01224202
## F 0.06682168 0.08220713 0.10275620 0.08849447 0.01224202 1.00000000statistical significance of the correlation matrix
colnames(cor_pvalues$p) <- cnames
rownames(cor_pvalues$p) <- cnames
print("and its statistical significance:")## [1] "and its statistical significance:"print(cor_pvalues$p)## A B C D E
## A 0.000000e+00 1.662902e-283 4.154338e-79 1.760060e-34 6.462174e-94
## B 1.662902e-283 0.000000e+00 0.000000e+00 3.566732e-226 6.802040e-04
## C 4.154338e-79 0.000000e+00 0.000000e+00 4.933289e-299 1.153638e-198
## D 1.760060e-34 3.566732e-226 4.933289e-299 0.000000e+00 0.000000e+00
## E 6.462174e-94 6.802040e-04 1.153638e-198 0.000000e+00 0.000000e+00
## F 7.410103e-09 1.417765e-45 5.378597e-38 1.399490e-54 0.000000e+00
## F
## A 7.410103e-09
## B 1.417765e-45
## C 5.378597e-38
## D 1.399490e-54
## E 0.000000e+00
## F 0.000000e+00acolor <- colorRampPalette(c("red", "white", "blue"))(20)
corrplot(
cor_matrix,
p.mat=cor_pvalues$p,
method="ellipse",
order="original",
col=acolor,
type='upper',
sig.level=0.0,
insig="p-value",
title='correlation matrix of original data'
)
Correlation summary of original data:
library(PerformanceAnalytics)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendchart.Correlation(dat1_noid, histogram=TRUE, pch=19)## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter
Correlation summary of detrended data:
library(PerformanceAnalytics)
chart.Correlation(dat1_noid, histogram=TRUE, pch=19)## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "method" is
## not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "method" is not a
## graphical parameter## Warning in plot.window(...): "method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "method" is not a graphical parameter## Warning in title(...): "method" is not a graphical parameter
plot covariance matrix of original data (pairwise)
cnames <- names(dat1_noid)
Ncnames <- length(cnames)
amatr <- matrix(unlist(dat1_noid), ncol=Ncnames, byrow=T)
colnames(amatr) <- cnames
cov_matrix <- cov(amatr)the covariance matrix of original columns set:
print("covariance matrix of original data:")## [1] "covariance matrix of original data:"print(cov_matrix)## A B C D E F
## A 5344.077 5082.750 4952.304 4846.738 4720.588 4654.519
## B 5082.750 5362.658 5149.259 4971.203 4812.842 4729.721
## C 4952.304 5149.259 5432.363 5141.714 4946.124 4834.244
## D 4846.738 4971.203 5141.714 5421.468 5071.708 4890.835
## E 4720.588 4812.842 4946.124 5071.708 5331.006 5071.083
## F 4654.519 4729.721 4834.244 4890.835 5071.083 5304.439p <- ggplot_heatmap(
inp=cov_matrix,
plot.title='Covariance matrix of original data'
)
print(p)
plot covariance matrix of detrended data (pairwise) Covariance is similar to variance. The latter applies to a single variable Covariance applies to two variables and is the sum of the product of each data-point’s deviation from the mean \[\sum_{i=1}^{N} \frac{(x_1-\bar(x_1))(x_2-\bar(x_2))}{N-1}\]
cnames <- names(dat1_detrended)
Ncnames <- length(cnames)
amatr <- matrix(unlist(dat1_detrended), ncol=Ncnames, byrow=T)
colnames(amatr) <- cnames
cov_matrix <- cov(amatr)the covariance matrix of detrended columns set:
print("covariance matrix of detrended data:")## [1] "covariance matrix of detrended data:"print(cov_matrix)## A B C D E F
## A 493.24296 -81.060207 -42.32725 -26.38814 -47.418538 -13.80354
## B -81.06021 548.650378 -111.09997 -72.66176 -8.280808 -35.65198
## C -42.32725 -111.099974 541.27997 -82.95156 -72.493255 -32.21031
## D -26.38814 -72.661759 -82.95156 496.51207 -112.577544 -37.23274
## E -47.41854 -8.280808 -72.49326 -112.57754 570.366279 -154.18776
## F -13.80354 -35.651978 -32.21031 -37.23274 -154.187763 608.91611p <- ggplot_heatmap(
inp=cov_matrix,
plot.title='Covariance matrix of detrended data'
)
print(p)
plot the original data
library(ggplot2)
library(reshape2)
adat = dat1_noid; adat[['id']] = dat1_id[['id']]
for(acol in names(dat1_noid)){
atitle <- paste0("original data, column ", acol)
cat("plot of ", atitle, "\n", sep='')
print(
ggplot(adat, aes_string(x='id',y=acol))+
geom_line()+
ggtitle(atitle)+
xlab('time')+
geom_smooth(method='gam')
)
}## plot of original data, column A
## plot of original data, column B
## plot of original data, column C
## plot of original data, column D
## plot of original data, column E
## plot of original data, column F
plot detrended data
adf <- data.frame(dat1_detrended); adf$id <- dat1_detrended_id[['id']]
for(acol in names(dat1_detrended)){
atitle <- paste0("detrended data, column ", acol)
cat("plot of ", atitle, "\n", sep='')
print(
ggplot(dat1, aes_string(x='id',y=acol))+
geom_line()+
ggtitle(atitle)+
xlab('time')+
geom_smooth(method='gam')
)
}## plot of detrended data, column A
## plot of detrended data, column B
## plot of detrended data, column C
## plot of detrended data, column D
## plot of detrended data, column E
## plot of detrended data, column F
discretise the data to five levels dat1_discrete <- discretise(dat1_noid, levels=5) plot discrete data
if( FALSE ){
for(acol in names(dat1_discrete)){
atitle <- paste0('Original data, discretised to 5 levels')
adf <- data.frame(data=dat1_discrete[[acol]]$discretised)
adf$id <- seq.int(nrow(adf))
p <- ggplot(
adf,
aes_string(x='id',y='data')
)+geom_line()+ggtitle(atitle)
print(p)
}
}discretise the detrended data dat1_detrended_discrete <- discretise(dat1_detrended, levels=5) plot detrended discrete data
if( FALSE ){
for(acol in names(dat1_detrended_discrete)){
atitle <- paste0('Detrended data, discretised to 5 levels')
adf <- data.frame(data=dat1_detrended_discrete[[acol]]$discretised); adf$id <- seq.int(nrow(adf))
p <- ggplot(
adf,
aes_string(x='id',y='data')
)+geom_line()+ggtitle(atitle)
print(p)
}
} # if( FALSE )find unique values
uniq_vals <- unique_values(inp=dat1_noid)
histo_dat1 <- histo(inp=dat1_noid,unique_vals=uniq_vals)
for(acol in names(dat1_noid)){
cat("found ", uniq_vals$num_unique_values[1,acol], " unique values in the original data, column ", acol, "\n", sep='')
barplot(
table(dat1_noid[[acol]]),
main=paste0('Original data, column ', acol, ' : ', uniq_vals$num_unique_values[1,acol], ' unique values.', sep=''),
ylab='unique levels'
)
}## found 261 unique values in the original data, column A
## found 65 unique values in the original data, column B
## found 38 unique values in the original data, column C
## found 656 unique values in the original data, column D
## found 958 unique values in the original data, column E
## found 63 unique values in the original data, column F
#print(histo_dat1)find unique values in the detrended data
uniq_vals_detrended <- unique_values(inp=dat1_detrended)
histo_dat1_detrended <- histo(inp=dat1_detrended,unique_vals=uniq_vals_detrended)
for(acol in names(dat1_detrended)){
cat("found ", uniq_vals_detrended$num_unique_values[1,acol], " unique values in the detrended data, column ", acol, "\n", sep='')
barplot(
table(dat1_detrended[[acol]]),
main=paste0('Detrended data, column ', acol, ' : ', uniq_vals_detrended$num_unique_values[1,acol], ' unique values.', sep=''),
ylab='unique levels'
)
}## found 76 unique values in the detrended data, column A
## found 30 unique values in the detrended data, column B
## found 35 unique values in the detrended data, column C
## found 45 unique values in the detrended data, column D
## found 941 unique values in the detrended data, column E
## found 175 unique values in the detrended data, column F
#print(histo_dat1_detrended)plot the histogram of each column
for(acol in names(dat1_noid)){
x <- dat1_noid[[acol]]
h <- hist(
x=x,
main=paste0("Histogram of original data, col ", acol, sep=''),
breaks=uniq_vals$histograms[[acol]]$breaks
)
#' and plot a normal density (if it was coming from that), from https://www.statmethods.net/graphs/density.html)
xfit<-seq(min(x),max(x),length=40)
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)
}
for(acol in names(dat1_detrended)){
x <- dat1_detrended[[acol]]
h <- hist(
x=x,
main=paste0("Histogram of detrended data, col ", acol, sep=''),
breaks=uniq_vals_detrended$histograms[[acol]]$breaks
)
#' and plot a normal density (if it was coming from that), from https://www.statmethods.net/graphs/density.html)
xfit<-seq(min(x),max(x),length=40)
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)
}
Density plots:
for(acol in names(dat1_noid)){
plot(
density(dat1_noid[[acol]]),
main=paste0("Density plot of original data, col ", acol, sep='')
)
}
comments on the distribution of the generating process of the original data: A : looks like a mixture of Gaussians (maybe 4) B : a lognormal/weibull distrbution C : a mixture of lognormal/weibull distributions or 1 lognormal/weibull and some normal distributions D : a mixture of normal distributions E : a lognormal/weibull (shape ~1) F : a mixture of lognormal/weibull and 1 normal(?)
for(acol in names(dat1_detrended)){
plot(
density(dat1_detrended[[acol]]),
main=paste0("Density plot of detrended data, col ", acol, sep='')
)
}
comments on the distribution of the generating process of the original data: A,B,C,F : mixture of normal distributions D,E : pure normal distributions with very light tails, see below on fitting a distribution on them Let’s also do a q-q plot (quantile-quantile plot) comparing with the Normal distribution. A q-q plot is a plot of the quantiles of test data coming e.g. from the Normal distribution against the corresponding quantiles of our data. ‘the data in the x% quantile’ means the first x% of the data sorted in ascending order. The more the plot coincides with the 45-degree line the stronger the assumption that our data derives from specific distribution is. Other distributions can be tested. note: flat regions indicate agreemenet, curved regions disagreement for example, flat in the middle but curved on the sides means that it probably comes from a Normal distribution but has heavy tails, more data at the tails than Normal.
library(ggpubr)
for(acol in names(dat1_noid)){
# qqnorm( y=dat1_noid[[acol]],
# main=paste0("Q-Q plot of original data, col ", acol, sep='')
# ); qqline(dat1_noid[[acol]]);
ggqqplot(
dat1_noid[[acol]],
main=paste0("Q-Q plot of original data, col ", acol, sep='')
)
}
for(acol in names(dat1_detrended)){
# qqnorm( y=dat1_detrended[[acol]],
# main=paste0("Q-Q plot of detrended data, col ", acol, sep='')
# ); qqline(dat1_detrended[[acol]])
ggqqplot(
dat1_detrended[[acol]],
main=paste0("Q-Q plot of original data, col ", acol, sep='')
)
}One can also try to fit a distribution to our data for example, in the detrended column D (which appears to have a pure normal distribution)
library(fitdistrplus)## Loading required package: MASS## Loading required package: survivalacol='D'
x <- dat1_detrended[[acol]]
afit <- fitdistr(x, 'normal')
print(afit)## mean sd
## -2.018302e-06 4.843524e-03
## ( 2.110433e-05) ( 1.492302e-05)plot(
density(x),
main=paste0('Density plot of detrended data, col ', acol, ' superimposed with fitted Normal dist.', sep='')
)
lines(
density(
rnorm(
n=50000,
mean=afit$estimate['mean'],
sd=afit$estimate['sd']
)
),
col='red',
pch=22
)
the q-q plot
qqplot(
x=x,
y=rlnorm(n=50000, afit$estimate),
main=paste0('Q-Q plot of original data, col ', acol, sep='')
)
In this case (column D), we can see the peak of the data is much higher than that which comes from the fitted distribution Here is column E (original) which appears as coming from lognormal distribution
acol='E'
x <- dat1_noid[[acol]]+0.02 #' must be positive!
afit <- fitdistr(x, 'lognormal')
print(afit)## meanlog sdlog
## 0.38862519 4.23427355
## (0.01844951) (0.01304578)and plot
plot(
density(x),
main=paste0('Density plot of original data, col ', acol, ' superimposed with fitted Normal dist.', sep='')
)
lines(density(rlnorm(n=50000, afit$estimate)), col='red', pch=22)
the q-q plot
qqplot(
x=x,
y=rlnorm(n=50000, afit$estimate),
main=paste0('Q-Q plot of original data, col ', acol, sep='')
)
which is not bad except the tail is too heavy! For those columns that I think derive from mixture processes we can use a mixture model estimator
library(mclust)## Package 'mclust' version 5.4
## Type 'citation("mclust")' for citing this R package in publications.acol='C'
x <- dat1_detrended[[acol]]
afit <- densityMclust(x)
plot(
afit,
what="density",
data=x,
breaks=50,
main=paste0('Fitting a mixture of Gaussian distributions to detrended data, col, ', acol, sep='')
)
and here is a q-q plot of original and fitted distribution notice that drawing samples from mclust fit is handled by function sample_from_mclust() in library lib/MIXTURES.R
qqplot(
x,
sample_from_mclust(mclustobj=afit, N=10000),
main=paste0('Q-Q plot of 10,000 samples from the fitted distr. and detrended data, col, ', acol, sep='')
)
Finally do a multi-variate fit using mclust. We used only 2 columns as it is better to visualise. All variables could have been used.
acol1='A'; acol2='B'
dat2 = cbind(
dat1_detrended[[acol1]],
dat1_detrended[[acol2]]
)
afit <- densityMclust(dat2)
plot(
afit,
what = "density",
type = "persp",
main=paste0('Multi-variate mixture fit for detrended data, cols ', acol1, ',', acol2, '.', sep='')
)
A mixture model for detrended column C with 6 Gaussians using package Rmixmod.
library(Rmixmod)## Loading required package: Rcpp## Rmixmod version 2.1.1 loaded
## R package of mixmodLib version 3.2.2
##
## Condition of use
## ----------------
## Copyright (C) MIXMOD Team - 2001-2013
##
## MIXMOD is publicly available under the GPL license (see www.gnu.org/copyleft/gpl.html)
## You can redistribute it and/or modify it under the terms of the GPL-3 license.
## Please understand that there may still be bugs and errors. Use it at your own risk.
## We take no responsibility for any errors or omissions in this package or for any misfortune that may befall you or others as a result of its use.
##
## Please report bugs at: http://www.mixmod.org/article.php3?id_article=23
##
## More information on : www.mixmod.orgacol='C'
nGaussians=6
x <- data.frame(data=(dat1_detrended[[acol]] - min(dat1_detrended[[acol]]) + 0.02))
afit <- mixmodCluster(x, nbCluster=nGaussians)
summary(afit)## **************************************************************
## * Number of samples = 52672
## * Problem dimension = 1
## **************************************************************
## * Number of cluster = 6
## * Model Type = Gaussian_pk_Lk_C
## * Criterion = BIC(-719436.1661)
## * Parameters = list by cluster
## * Cluster 1 :
## Proportion = 0.1782
## Means = 22.0200
## Variances = 0.0000
## * Cluster 2 :
## Proportion = 0.1782
## Means = 22.0200
## Variances = 0.0000
## * Cluster 3 :
## Proportion = 0.1556
## Means = 18.5940
## Variances = 3.7390
## * Cluster 4 :
## Proportion = 0.1704
## Means = 20.6881
## Variances = 0.2275
## * Cluster 5 :
## Proportion = 0.1633
## Means = 23.3716
## Variances = 0.2362
## * Cluster 6 :
## Proportion = 0.1543
## Means = 25.5165
## Variances = 3.8131
## * Log-likelihood = 359810.4937
## **************************************************************histCluster(
x=afit["bestResult"],
data=x,
breaks=100,
# CHECK THIS:
main=c(paste0('Histogram of detrended data, col ', acol, ' fitted with ', nGaussians, ' Gaussians.', sep=''))
)
Here we estimate a mixture of different distribution families For example Gaussian and Gamma for original data, column C
library(flexmix)## Loading required package: latticeadat <- list()
acol='C'
adat[[acol]] = dat1_noid[[acol]]-min(dat1_noid[[acol]])+10we will a mix of 3 Gaussians and 1 Gamma
mymodels <- list(
FLXMRglm(family = "gaussian"),
FLXMRglm(family = "gaussian"),
FLXMRglm(family = "gaussian"),
FLXMRglm(family = "gaussian")
)use the data from column C (original) for k=6 components
st<- system.time(
afit <- flexmix(
C ~ 1,
data=adat,
k=6,
model=mymodels
)
)
cat("done, model estimated in ", st[2], " seconds.\n", sep='')## done, model estimated in 2.535 seconds.print(parameters(afit))## [[1]]
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## coef.(Intercept) 27.680767 10.2295131 15.9485114 20.347395 13.1970149
## sigma 3.602226 0.4205237 0.8338333 1.704873 0.3977475
## Comp.6
## coef.(Intercept) 1.200000e+01
## sigma 1.198945e-12
##
## [[2]]
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## coef.(Intercept) 27.680767 10.2295131 15.9485114 20.347395 13.1970149
## sigma 3.602226 0.4205237 0.8338333 1.704873 0.3977475
## Comp.6
## coef.(Intercept) 1.200000e+01
## sigma 1.198945e-12
##
## [[3]]
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## coef.(Intercept) 27.680767 10.2295131 15.9485114 20.347395 13.1970149
## sigma 3.602226 0.4205237 0.8338333 1.704873 0.3977475
## Comp.6
## coef.(Intercept) 1.200000e+01
## sigma 1.198945e-12
##
## [[4]]
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## coef.(Intercept) 27.680767 10.2295131 15.9485114 20.347395 13.1970149
## sigma 3.602226 0.4205237 0.8338333 1.704873 0.3977475
## Comp.6
## coef.(Intercept) 1.200000e+01
## sigma 1.198945e-12obtain additional information on the estimated model using refit()
arfit <- refit(afit)## Warning in sqrt(diag(z@vcov)[indices]): NaNs produced
## Warning in sqrt(diag(z@vcov)[indices]): NaNs produced
## Warning in sqrt(diag(z@vcov)[indices]): NaNs produced
## Warning in sqrt(diag(z@vcov)[indices]): NaNs producedplot(
arfit,
bycluster=F,
main=paste0('Mixture Model Components for orig. data, col ', acol, sep='')
)
compare densities, sample from our estimated model using rflexmix() :
sam <- rflexmix(afit)
plot(
density(adat[[1]]),
main=paste0('Fitting a mixture of different distributions to detr. data, col ', acol, sep='')
)
lines(density(sam$y[[1]]), col='red', pch=22)
not a bad fit at all!